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This article is aimed at extending the framework of optimal control techniques to the situation 
where the control field values are restricted to a finite set. We propose a generalization of the 
standard GRAPE algorithm suited to this constraint. We test the validity and the efficiency of 
this approach for the inversion of an inhomogeneous ensemble of spin systems with different offset 
frequencies. It is shown that a remarkable efficiency can be achieved even for a very limited number 
of discrete values. Some applications in Nuclear Magnetic Resonance are discussed. 
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I. INTRODUCTION a schematic description). 


The design of control sequences accounting for exper¬ 
imental constraints is a central task in a variety of do¬ 
mains in quantum dynamics extending from photochem¬ 
istry, Nuclear Magnetic Resonance (NMR) and quantum 
information science M- Nowadays, Optimal Control 
Theory (OCT) reveals to be a highly efficient and ver¬ 
satile tool to bring answers to the different issues raised 
by the experimental setups For the past few 

years, there has been an intense theoretical activity in 
developing new optimal control procedures able to build 
high quality control fields in presence of some experimen¬ 
tal imperfections and constrain ts 111 , [l3l l2dl - !25[ . These in¬ 
clude spectral constraints |26 h 29| . amplitude and phase 
transients [30| . non-linear interactions between the sys¬ 
tem and the control field [3ll - l33j . robustness against ex¬ 
perimental uncertainties and errors [2C| [HJ . To date, the 
majority of studies has assumed that the amplitude and 
phase of the control field can vary continuously. For ex¬ 
ample, in modern NMR spectrometers |7], the amplitude 
and phase of the control pulses can be defined with high 
resolution, allowing for a virtually continuous variation 
of these parameters [HI, [36|. More generally, this is pos¬ 
sible in experimental sett ing s, where arbitrary waveform 
generators are available [301 1371 - 1391 ]. However, in many 
cases the available hardware only allows to switch be¬ 
tween a discrete set of pulse phases [4(], EJ . This set of 
phases can be chosen before the experiment, but cannot 
be altered during the experiment. Hence the control is 
quantized and restricted to a fixed finite number of val¬ 
ues, which can nevertheless be optimized (See Fig. |T] for 
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FIG. 1. Schematic representation of the mapping vector be¬ 
tween the set of M- phases and the set of N- time steps. The 
mapping is depicted by the dashed lines. 

This paper is aimed at tackling this problem by propos¬ 
ing two fundamentally different approaches. The first 
method is based on an extension of a standard optimal 
iterative procedure, namely GRAPE (HI El. The sec¬ 
ond option can be viewed as a brute force strategy in 
the sense that a standard optimization method is used 
in a first step, producing a continuous control field. A 
Lloyd procedure, analogous in its spirit to a mean square 
method, is then applied to quantize the field [43|. We test 
the validity and the efficiency of the two approaches on a 
benchmark control problem, the simultaneous control of 
an ensemble of uncoupled spin 1/2 particles with different 
offset terms 13,1H| . Extensive numerical simulations 
reveal that an efficient control can be achieved even for a 
small number of quantized values of the control field. We 
analyze the relative efficiency of the two algorithms as a 
function of the number of allowed values for the control 
field. Finally, note that the discretization of the control 
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field has been recently used in a series of paper to ac¬ 
celerate the numerical integration of the time-dependent 
Schrodinger equation |4dH50l |. This idea based on the 
precomputation of elementary propagators was applied 
to quantum optimal control problem and would be also 
useful in our case. 

The remainder of the paper is organized as follows. In 
Sec. nn we outline the principles of the new optimization 
procedures, paying special attention to the flexibility and 
to the applicability of the methods. Section [TIT] is dedi¬ 
cated to the presentation of the numerical results. The 
efficiency of the quantized control field is compared to its 
continuous counterpart. The advantages of the optimiza¬ 
tion algorithm directly accounting for the constraint and 
of the brute force strategy are also discussed. Conclusion 
and prospective views are given in Sec. IIVI 

II. THEORY 

The goal of this section is to propose an optimal con¬ 
trol algorithm suited to quantum systems where the con¬ 
trol field is restricted to a finite set A4 of M discrete 
values M. = ,%}, where M is fixed and the 

values Vk are optimized. To simplify the discussion, we 
consider here the case of mixed quantum systems, which 
also applies directly to pure states and can be straight¬ 
forwardly extended to the control of unitary operators in 
quantum computing. Starting from an initial state po, 
optimal control is invoked in order to maximize the pro¬ 
jection onto a target state p/, which is measured by the 
following figure of merit: 


be straightforwardly adapted to phase-modulated pulses 
with a constant amplitude [58], as used in Sec. OH] 

Let U ( t ) be the evolution operator associated with the 
Hamiltonian H(t), such that p(t) = U[t)p$U\t). Fol¬ 
lowing the description of a standard GRAPE algorithm 
[25l ] , we assume that the field u is discretized in time with 
a time step At = tf/N , where N is the number of val¬ 
ues of the field. We denote by Uj : with j = 1,2, ■ • • ,7V 
the value of u in the interval [(j — l)At, jAt\. In this 
approximation, note that the control field u is now de¬ 
scribed by a set of TV reals (tti,U2 , ••• ,ujv)- The time 
evolution can be computed by the TV block propagators 
(Ui, U 2 , ■ ■ ■ ,U N ), where 

JJj = e ~ iAt ( H o+ u j H i)' ( 4 ) 

The optimal control problem is described by a figure of 
merit $ to be maximized. Using the different evolution 
operators Uj. >I> can be written as follows: 

$ = Tr(U N U N -i ■ ■ ■ U 1P0 UM • • • Ulfpj). (5) 

The optimization procedure is based on the derivation 
of the gradient of the figure of merit with respect to the 
different values of the control field: 

H = Ty(U n U N -i •••§!••• IhpoUM • • • uUf) 

+ Tr (U N U N -! ■ ■ ■ U 1P0 U\UI Uy f ). 

( 6 ) 

We assume that the different propagators are approxi¬ 
mated by a first order split-operator, which is valid up to 
the order 2 in At: 

77. ~ 


$ = Tr [p\p(t f )\, (1) 

where tf is the control duration and p(t) the state of the 
system at time t. 


A. The GRAPE algorithm 


In this paper, we consider specifically the GRAPE al¬ 
gorithm [25|, but the same construction of the discrete 
version can be used for other al gor ithms such as the 


monotonic or Krotov ones 1321. 51 15 If (see also the gen¬ 


eral analysis of such methods [55], [56( ). 

The time evolution of the mixed-state p(t ) is ruled by 
the Liouville-von Neumann equation: 


ip(t) = [H(t),p(t)\, (2) 

where units such that H = 1 have been chosen. The 
Hamiltonian H(t) of the system can be written as: 

H{t) = H 0 + u(t)Hi, (3) 

where Hq is the field-free Hamiltonian operator and 
Hi the interaction part. A general presentation of 
the algorithm is proposed here, but the procedure can 


This hypothesis simplifies the computation of the deriva¬ 
tive of Uj with respect to Uj, which can be written as: 

r)JI- 

d^^Uji-iAtHi). ( 8 ) 

Note that a different approximation of the propagator 
would give a different derivative. A straightforward com¬ 
putation using Eq. ([5]) then leads to the gradient of the 
figure of merit: 

7 ^- - —iAtTi(Xj [H-[ , pj \), (9) 

where the states pj and Aare defined by: 

Pj = Uj-i ■ ■ • u lPo ul ■ • • 

*j = u]u} +1 .~uy f u N -u j . 

In a first-order GRAPE algorithm (see [57I ] for a recent 
generalization to the second order), the field u is updated 
at each step by the formula: 

u ‘^ Ui + e a^' <10) 

with e a small real parameter, which is chosen from a 
line search method to ensure the increase of the figure of 
merit $. 
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B. A discrete-pulse version of GRAPE 


In the continuous version of the GRAPE algorithm, 
the control field u can take any real value, i.e. Uj £ R. 
We consider in this section that u is restricted to a finite 
set At of M values. We introduce the mapping vector 
from the set {1, 2, • • • ,N} to {1, 2, • • ■ M} which as¬ 
sociates with the different values of u, a value in Ai: 
u j = v ll{j) ■ This mapping makes a correspondence be¬ 
tween the set of time steps and the set of dicrete values. 
A schematic illustration of this process is given in Fig. |Tj 
This optimization procedure has a non trivial charac¬ 
ter in the sense that both the discrete values {u m } and 
the mapping can be adjusted to maximize the figure 
of merit •!>. In the algorithm proposed in this paper, 
each iteration is decomposed into two sub-steps. In the 
first sub-step, the mapping is fixed and the {v m }~ val¬ 
ues are optimized through a modified version of GRAPE 
which can be described as follows. The functional $ can 
be derived with respect to v m : 


<9<f> 


dv 


m 


-iA<Tr(Aj[ffi,pj]). 


( 11 ) 


Note that the derivative of the quality factor with re¬ 
spect to the discrete value v m is simply the sum over the 
derivatives with respect to Uj which maps to v m . The 
control v m is then updated at each step of the algorithm 
through the formula: 

Vm 1 t Vm + £— . (12) 

The roles are reversed in the second stage which aims 
at optimizing p, while the discrete values {u m } are not 
changed. Since the total number N of values of the field 
can be very large, it would be time consuming to find the 
global optimal mapping . Instead, we adopt a more 
direct approach which allows to improve the final result 
without no guarantee to attain its upper bound. The 
efficiency of this procedure will be shown numerically and 
discussed in Sec. m 

For each time step j £ {1, 2, • • ■ , N} taken in increas¬ 
ing order, we test the M possible values of the control 
field Uj = vi, l = 1, 2, • • • , M, by computing the corre¬ 
sponding figure of merit $>(uj = vi). Then we define the 
new discrete phase as being the one which maximizes the 
quality factor: 

p(j) = k = arg max <F (rife) 

A proper use of the adjoint state allows us to obtain a 
computational cost of the algorithm with depends lin¬ 
early on M. 

Since our method is a two-step procedure, it cannot be 
easily extended to the second order (the second order is 
related here to the gradient and not to the order of accu¬ 
racy of the propagator). Note also that this approach is 
compatible with a toolkit method [461450^ . At each step 


of the algorithm, the propagators associated with the set 
of M- discrete values can be precomputed. These prop¬ 
agators can be used to update the mapping at iteration 
k — 1 and the same set of propagators allows us to derive 
the gradient at step k of the algorithm. 


C. Quantization and Lloyd’s algorithm 


We propose in this paragraph a second strategy based 
on Lloyd’s algorithm. The idea consists first in using a 
standard GRAPE algorithm to build a continuous control 
field. In our numerical example, we will be interested in 
a case where the phase of the field is optimized, while its 
amplitude is constant. Nevertheless, it would be straight¬ 
forward to extend this procedure to a general situation 
with no constraint on the control field. Lloyd’s approach 
is then applied in a second step to quantize this field. 
Lloyd’s algorithm, which is explained in detail below (we 
refer the reader to [43| for additional information), can be 
viewed as a mean-square approximation procedure which 
allows us to approach a continuous function by a discrete 
one. In contrast to the discrete version of GRAPE pre¬ 
sented in Sec. IIIB1 no information about the dynamics is 
used for the quantization. In other words, the computa¬ 
tion of the continuous control field is sufficient to design 
its discrete counterpart. However, since only a geomet¬ 
ric (and not a dynamical) information is used, there is a 
priori no guarantee about the efficiency of the quantized 
control. 

To simplify the description of the algorithm, we con¬ 
sider a set (iti)i=i,...,jv of values in the interval [0, 27 t[ (see 
Sec. (mi for details), which correspond to the pulse phases 
derived from the continuous GRAPE algorithm. The an¬ 
gles (ut) are defined modulo 27 t and we identify 0 and 27r 
by periodicity. The discrete control field values derived 
from Lloyd’s algorithm will be denoted We 

also introduce two sets of M reals belonging to [0, 27r[, 
B (k) = (Bj k )) j=1 "" iM and = (Y / k) ) j= sorted 
in increasing order, k being the iteration step of the algo¬ 
rithm. The initial set is defined such that the values 
Bj °' 1 are equally distributed in increasing order in the in¬ 
terval [0, 27 t[. Note that the 27 t- periodicity of the control 
field is taken into account in the different relations used 
in the algorithm, even if this point is not explicitly men¬ 
tioned below in order to clarify the presentation of the 
computation of the discrete field. The algorithm can be 
described as follows. 

Algorithm: Given the initial set B^°\ the control 
phases ('Ui);=i,... l Ar, and the iteration parameters k = 0, 
e > 0, Jo = 0, do: 

1. k = k + 1. 


2. For j = 1, • • • , M, do 
’■( k ) _ Mean^U, 


Y i 


TD (k—1) n(fc —1) 

’ ° 3 +1 


3. For * = !,■•• , N, evaluate 
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dt = miiij 


|Ui — Yj k ' i mod 2ir 

i di 


4. if \J k - Jk- 1 | > e, do 

For j = 1, • ■ • , M, B {k) = (Y} k) + YjU)/2 
and go to 1 
else go to 5 


5. For i = 1, • • • ,7V, compute uji = Proj(uj, F^) 

Note that the function ” Mean” stands for the mean value 
of a set of numbers belonging to a given interval. The 
function ”Proj” denotes the projection of each iq onto 
the discrete set of values Y W. More precisely, we define 
the projection as the value uj-i £ y( fe ) which minimizes 
the distance from tq to Y^ k \ It is then straightforward 
to define the mapping and the set A4 from the values 
such that v = u>». A schematic descrip¬ 
tion of this algorithm is displayed in Fig. [2] 


III. APPLICATION TO THE CONTROL OF 
SPIN SYSTEMS 

This section is dedicated to the application of the dis¬ 
crete version of the GR APE algorithm for controlling the 
dynamics of spin systems driven by radio-frequency mag¬ 
netic fields. The different numerical values are chosen so 
as to reproduce the typical features of a spin sample in 
liquid state NMR spectroscopy with Bq- inhomogeneities 
0. We test the validity and the efficiency of the two 
approaches developed above for the inversion of an en¬ 
semble of inhomogeneous uncoupled spins with different 
resonant offset frequencies [44 ]. 


A. The model system 

In a given rotating frame, the Hamiltonian of each 
isochromat, i.e. each sub-system with a different reso¬ 
nance frequency w, is H u = ua z + u) x cr x + u} y a y . For 
a matter of convenience and to simplify the numerical 
implementation of this problem, we move to the Bloch 
picture 0: 


= - OjMy + lu v M z (13) 

M^ ] =loM x - u x M x 
=UJ X My - 0J y M X , 



where the Bloch vector of the isochromat M ^ = 
(Mx^, My ^, Mj w )) represents the state of the system, 
which can be readily related to the corresponding den¬ 
sity matrix through the identification = Tr (p^a), 
with i = {x, y, z} and cr, the Pauli matrices. On the 
right hand side of Eq. C3D, uj is the offset frequency 
term which belongs to the interval [— u) max , uj m ax]- The 
two components of the control field along the x- and y- 


FIG. 2. (Color online) Schematic representation of the differ¬ 
ent steps of Lloyd’s algorithm in the case M = 4. The dots 
indicate the position of the control field values {ui}. The two 
sets of reals B ^ an d y( fc ) are represented by solid and dashed 
lines, respectively. At iteration k = 1, the top panel shows 
the way to compute the 1 * values according to the step 2 
of the algorithm. Different colors have been used to help the 
understanding of the procedure. Starting from this new set 
Y (1) , the set B (1) is defined in the middle panel by using step 
4 of the algorithm. The bottom panel depicts a new step 2 
for the next iteration of the algorithm. At this stage, note 
that the new Y- values have not been computed. 
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directions are denoted co x and co Vl respectively. Starting 
—y —y 

from M(0)=M z , the goal of the control is to reach the 
south pole of the Bloch sphere = — A% z for any spin 
of the ensemble. The quality factor or the figure of merit 
to maximize can be written as follows: 

$ =-V (14) 

n °ff , 

J J 1=1 

where n 0 ff is the total number of uncoupled spins and 
tf the total control time. The offset terms 04 are chosen 
equally spaced in [—w max , w maa: ]. Note that this control 
problem has been extensively investigated in the litera¬ 
ture in the standard situation where the field is continu¬ 
ous. We refer the interested reader to key publications on 
this subject f20j. Following Ref. [44j, the numerical val¬ 
ues are taken to be u; max /(27r) = 10 kHz, n 0 ff = 200 and 
tf = 0.18 ms. The control field at = (ui x ,uj y ) is assumed 
of fixed control amplitude uiq and can be expressed as 

at = cao [cos(0(f))"at + cos(0(f))l/], (15) 

where 9 £ [0,27r[ is the phase to optimize. The maximum 
pulse amplitude o;o/(27r) is chosen constant and equal to 
10 kHz, the same value as oj mQX . The time digitization is 
taken as Af=0.5 ps. We have checked that the qualitative 
conclusions of this paper do not depend on a specific 
choice of the used constants. 


the control field, but a significant number of examples 
converges towards a quality factor close to 1 , even for 
M = 4. For instance, Fig. [3] shows that a large percent¬ 
age of the optimizations is close to the maximum quality 
factors of 0.992, 0.993, 0.994 and 0.9965 for the cases of 
M = 4, M = 8 , M = 12 and M = 16, respectively. As 
could be expected, we observe that the higher the num¬ 
ber of discrete values M , the narrower the distribution is. 
Note also that the global shape of the histogram distri¬ 
bution in this discrete setting is very similar to standard 
distributions that can be observed in the continuous case 
[dll . l45l | . Another standard choice for the initial values 
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B. Numerical results 

We are now in a position to check the efficiency of the 
discrete GRAPE algorithm through the comparison with 
the continuous optimal solution. In all the numerical 
simulations with the continuous version, we choose for 
the initial control phase the following simple form: 

e(t) =f( J,- 1 ) 2 ' (16) 

This parabolic behavior leads to a linear evolution of the 
corresponding frequency. This class of control fields is 
known to be robust to experimental imperfections in the 
adiabatic limit (5f|. With this initial guess, the GRAPE 
algorithm converges to the target state with an accuracy 
better than 0.9982 and with an optimal control field very 
close to the solution proposed in Ref. 144] , see also Fig. 
[ 6 ] for a plot of this optimal field. 

In the discrete case, we recall that M is the number of 
possible discrete values that can be taken by the control 
phase 9. Figure [3] shows the histogram distribution of 
the quality factor $ for the cases M = 4, 8 , 12 and 16. 
In each situation, the initial discrete values are randomly 
chosen in the interval [ 0 , 27r[ to generate 100 possible re¬ 
alizations. The initial mapping p is also generated ran¬ 
domly. We observe that the final quality factor depends 
on the initially chosen (random) discrete set of values of 


FIG. 3. Histogram distribution of the final quality factor 
4> obtained with the discrete GRAPE method, for the cases 
M = 4 ((a)- red bars), M = 8 ((b)- blue bars), M = 12 ((c)- 
green bars) and M = 16 ((d)- black bars). In each optimal 
computation, the initial discrete values are randomly chosen 
in the interval [0, 2zr[ to generate 100 possible realizations. 

of the control field is an equidistant repartition in the 
interval [0,2-7 t[. Here, we consider an initial zero control 
field and we determine the initial mapping by a forward 
propagation where the M possible values of the field are 
tested for each time step. Note that this choice is not 
crucial for the final efficiency of the algorithm, but it 
allows in one shot to achieve a very high performance. 
The corresponding figures of merit are depicted in Fig. |4] 
as a function of the number of iterations. Here again, 
we observe the good behavior of the algorithm since a 
performance higher than 0.99 is achieved for M > 8 . 

As could be expected, the larger the number of discrete 
values M is, the higher the efficiency. This point is clearly 
shown in Fig. [5] where the value of the final quality factor 
<F is plotted for the two specific choices of the initial set 
of discrete values as a function of the number of discrete 
values M. 

For the case of a random choice of the initial set, we 
observe that the averaging over 100 realizations of the 
optimized quality factor seems less efficient than the per¬ 
formance achieved when the initial discrete values are 
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FIG. 4. (Color online) Plot of the quality factor <3? obtained 
with the discrete GRAPE method as a function of the number 
of iterations. For this example, the initial values of the phase 
are uniformly distributed in the interval [0,27r[. 


uniformly distributed in the interval [0, 2n[. However, 
the maximum of the hundred realizations have a very 
high performance, which is better than the one with a 
uniform distribution, even for M = 4. This point is illus¬ 
trated in Fig. 0 where the highest quality factor among 
the 100 random realizations is plotted as a function of 
M. This observation illustrates a standard optimization 
problem with a gradient procedure, which only gives a lo¬ 
cal information and cannot avoid traps and local maxima 
in the control landscapes |.2j. 

The complexity of the optimal solutions designed by 
the algorithm is illustrated in Fig. [G] for M = 4, 6 and 
8 . For these particular cases, we get an efficiency of the 
order of 0.99 for M = 4, 0.9925 for M = 6 and 0.9931 
for M = 8. 

In a second series of simulations, we explore the ef¬ 
ficiency of the Lloyd’s algorithm for the same control 
problem. As explained in Sec. Ill Cl we first apply a 
standard continuous version of GRAPE. Figure [6] shows 
the different discrete fields designed by this procedure, 
which remain very close to the continuous one. Note the 
difference with respect to the solution determined from 
GRAPE in Fig. [G] The representation used in Fig. 0 
helps the visualization of the two different results. 

The numerical findings of Fig.[8]confirm the superiority 
of the discrete version of GR APE over Lloyd’s algorithm 
when M < 6. For larger values of M, we observe that 
the Lloyd’s algorithm becomes more efficient. This point 
can be interpreted in the light of the local procedure used 
here to design the mapping. Optimizing such a mapping 
is intrinsically a global combinatorial optimization prob¬ 
lem which cannot be easily replaced by a local search 
method. Another choice could be to combine a GRAPE 
approach, together with a heuristic global procedure such 
as genetic algorithms, which would be used to optimize 


FIG. 5. (Color online) Plot of the final quality factor <f> for 
two specific initial sets of discrete values as a function of the 
number of discrete values M. In the case of a random initial 
set (shown by the square blue solid line), the quality factor 
has been obtained by averaging over the 100 realizations given 
in Fig. [3] The red solid line with circles represents the value 
of the final figure of merit $ when the initial discrete values 
of the control field are uniformly distributed in the interval 
[0, 27 t[. The black curve with diamonds displays the highest 
quality factor among the 100 random realizations. 
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FIG. 6. (Color online) Examples of optimal discrete control 
phases obtained with the discrete GRAPE method (left col¬ 
umn) and with the Lloyd’s algorithm (right column) for the 
cases M = 4 (a-d), M = 6 (b-e) and M = 8 (c-f). For 
GRAPE, the initial discrete values are randomly chosen in 
the interval [0, 2-7r[. The continuous optimal field is plotted in 
solid line (red) in the right panels. 
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the mapping Q • The global nature of the search method 
seems interesting even if the computational cost associ¬ 
ated with this approach may be prohibitive. 


As could be expected, the figure of merit for the Lloyd 
procedure tends to the continuous optimal result as the 
number of discrete values goes to infinity. This is not the 
case for the discrete version of GRAPE when random 
or uniform initial phases are used. In order to improve 
the efficiency of the quantization procedure, we have also 
used the values derived by Lloyd’s algorithm as initial 
guess for GRAPE. However, we observe, for this example, 
that the combination of the two algorithms does not lead 
to a significantly better figure of merit. Again this is 
probably a direct consequence of the approach we choose 
to optimize the mapping. The comparison of the relative 
performance of the two algorithms should be tested on 
other control problems to confirm the conclusions of this 
paper. 


© 
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FIG. 7. Representation on a circle of the discrete phases com¬ 
puted in the example of Fig. [fj] for the GRAPE (left column) 
and the Lloyd (right column) algorithms with M = 4 (top, 
a-d), M = 6 (middle, b-e) and M = 8 (bottom, c-f). 


IV. CONCLUSION 

This paper has focused on the application of new 
optimization procedures for studying the simultaneous 
control of an ensemble of uncoupled spins with different 
offsets. A basic feature of the methods under considera¬ 
tion is that they allow us to account for a quantization of 
the control, which can only take a fixed finite number of 
values. Numerical tests have demonstrated the efficiency 
of the proposed procedures, showing that a high quality 
control can be achieved even with a very small number 



FIG. 8. (Color online) Comparison of the figures of merit 
achieved with the discrete version of GRAPE and Lloyd’s al¬ 
gorithm for different numbers of discrete values. In the case of 
GRAPE, the best result over the 100 realizations of Fig.[3]has 
been used. The horizontal solid line indicates the efficiency 
of the continuous version of GRAPE. 

of quantized values. This conclusion is crucial in some 
electronically controlled nanodevices or in Nuclear 
Magnetic Resonance for which such constraints have to 
be taken into account. Experiments are in progress in 
the field of Electron Paramagnetic Resonance in which 
such constraints are imposed by the hardware [60j. 
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